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We present both spin-independent and -dependent parts of a central interquark potential for 
charmonium states, which is calculated in 2-1-1 flavor dynamical lattice QCD using the PACS-CS 
gauge configurations with a lattice cutoff of ~ 2.2 GeV. Our simulations are performed with a 
relativistic heavy-quark action for the charm quark at the lightest pion mass, — 156(7) MeV, in 
a spatial volume of (3 fm)^. We observe that the spin-independent charmonium potential obtained 
from lattice QCD with almost physical quark masses is quite similar to the Cornell potential used 
in nonrelativistic potential models. The spin-spin potential, which is calculated in full lattice QCD 
for the first time, properly exhibits a finite-range repulsive interaction. Its r-dependence is different 
from the Fermi-Breit type potential, which is widely adopted in quark potential models. 



Recently, many of the newly discovered charmonium- 
like mesons have been announced by B-factories at KEK 
and SLAC, which are primarily devoted to the physics 
of CP violation. Such states named as "XYZ" mesons 
could not be simply explained by a constituent quark 
description in quark potential models [T]. Thus, the 
charmoniumlike XYZ mesons are expected to be good 
candidates for nonstandard quarkonium mesons such as 
hadronic molecular states, diquark-antidiquark bound 
states (tetraquark states) or hybrid mesons How- 
ever, it seems to be still too early to judge whether we 
may discard the constituent quark description for such 
higher-mass charmonium states. 

The interquark potential in quark potential models 
is based on the phenomenology of confining quark in- 
teractions: the so-called Cornell potential [3] together 
with spin-spin, tensor and spin-orbit terms appeared as 
leading spin-dependent corrections in powers of the in- 
verse of the heavy-quark mass mg [U Although 
the Cornell potential was qualitatively justified by the 
static heavy-quark potential obtained from Wilson loops 
in lattice QCD [5], the functional forms of the spin- 
dependent terms in potential models are basically deter- 
mined by perturbative one-gluon exchange as the Fermi- 
Breit type potential [4 . Thus, properties of higher-mass 
charmonium states predicated in potential models may 
suffer from large uncertainties, since the phenomcnologi- 
cal spin-dependent potentials would have validity only at 
short distances and also in the vicinity of the heavy-quark 
mass limit. 

In this sense, the reliable interquark potential directly 
derived from first principles QCD is desired at the charm 
quark mass. Indeed, the static potential between in- 
finitely heavy-quark and antiquark, which is obtained 
from Wilson loops, have been precisely calculated in 
lattice QCD simulations [S]. The relativistic correc- 
tions to the static potential are classified in powers of 
l/mg within a framework called potential nonrelativis- 
tic QCD [7j. The leading and next-to-leading order cor- 



rections have been successfully calculated in quenched 
lattice QCD with high accuracy by using multilevel algo- 
rithm |8]. However, it is rather difficult to calculate the 
proper charmonium potential in lattice QCD within the 
Wilson loop formalism. It is obvious that the inverse of 
the charm quark mass is far outside the validity region 
of the l/mg expansion. Indeed, a spin-spin potential de- 
termined at ©(l/mg) [Hj, which exhibit an attractive in- 
teraction for the higher spin states, seems to yield wrong 
mass ordering among hyperfine multiplets. In addition, 
practically, the multilevel algorithm employed in Ref. [5] 
is not easy to be implemented in dynamical lattice QCD 
simulations. 

Under these circumstances, in our previous work fS], 
we have proposed a novel approach, where the interquark 
potential at finite quark mass can be accurately deter- 
mined from the equal-time and Coulomb gauge Bethe- 
Salpeter (BS) amplitude through an effective Schrodinger 
equation (See also Ref. (TU]). The BS amplitude method 
is originally applied for the hadron-hadron interaction 
jll| . As demonstrated in Ref. [9], the spin- independent 
part of an interquark potential calculated in the new 
method smoothly approaches the static potential given 
by Wilson loops in the infinitely heavy-quark limit. 
The new approach enables us to determine both spin- 
independent and spin-dependent interquark potentials at 
the charm quark mass, which potentially account for all 
orders of l/mg corrections. Furthermore, there is no 
restriction to dynamical calculation within this method. 

In this paper, we present results of the spin- 
independent central and spin-spin potentials between the 
quark (Q) and antiquark (Q) at the vicinity of the phys- 
ical charm quark mass from the BS amplitude with 2-1-1 
fiavor PACS-CS gauge configurations [12], where the 
simulated pion mass is closest to the physical point as 
Mjr — 156(7) MeV. As for a treatment of heavy quarks, 
we adopt a relativistic heavy quark (RHQ) action, which 
can remove large discretization errors introduced by large 
quark mass [13'. 
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TABLE I: Summary of RHQ parameters and results of S- and P-wave charmonium masses. 
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cb 
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B -r6'-wavc 


Tt T^'-wavc 
J" hyp 


r). (0-+) 


J/^ (1-+) 


Xco (0++) 


Xcl (1 + +) 


/ic (1+-) 


[GeV] 


3.0638(9) 


0.1124(9) 


2.9794(5) 


3.0919(10) 


3.3865(58) 


3.4781(62) 


3.4995(62) 



Let us briefly review the new method utihzed here to 
calculate the interquark potential with the finite quark 
mass. As the first step, we consider the following equal- 
time QQ BS amplitude in the Coulomb gauge for the 
quarkonium states [TH [T3] : 

(/.r(r) = 5](0|Q(x)rO(x + r)|QQ; J^^), (1) 

X 

where r is the relative coordinate of two quarks at certain 
time slice and F is any of the 16 Dirac 7 matrices. A sum- 
mation over spatial coordinates x projects onto the zero 
total momentum. The r-dependent amplitude, (/)r(r), is 
here called BS wave function. 

In lattice simulations, the BS wave function can be ex- 
tracted from the following four-point correlation function 

^ (0|Q(x, i)rg(x + r, t) (g(x', ts)rQ(y', ts))^ |0) 

x,x',y' 

^ E E A.(0|Q(x)rQ(x 4- r) In) e-*^" 

>Ao'/>r(r)e-<(*-*=), (2) 



X n 



at the large Euclidean time from source location (tg)- 
Here both quark and antiquark fields at are sepa- 
rately averaged in space as wall sources. The constant 
amplitude is a matrix element defined as = 

Ex'.y'("l (Q(x')ri3(y'))^ |0)- ^1 denotes a mass of the 
n-th quarkonium state \n) in a given J^*^ channel. For 
instance, when F is chosen to be 75 for the pseudoscalar 
(PS) channel {J^^ = 0~+) and 7^ for the vector (V) 
channel {J^^ = 1 ) in the charm sector, Mq^ and 
correspond to the rest masses of the 7/c and J /ip ground 
states, which can be read off from the asymptotic large- 
time behavior of the correlation functions. 

The BS wave function satisfies an effective Schrodinger 
equation with a nonlocal and energy-independent in- 
terquark potential U [HI [16] 



V2 
2/i 



0r(r) + j dr'[/(r,r')^r(r') = ^r-^rW, 



(3) 



where the reduced mass ji of the QQ system is given 
by mQ/2 with the quark kinetic mass mg. The energy 
eigenvalue E-p of the stationary Schrodinger equation is 
supposed to be Afp — 2mQ [17^ . If the relative quark ve- 
locity V — iV/mgl is small as v <^ 1, the nonlocal poten- 
tial U can generally expand in terms of the velocity v as 



U{v\v) = {y(r) + l/s(r-)SQ-SQ + yT(r)5i2 + VLs(r)L-S + 
0{v^))5{v' - r), where = (Sq • f)(SQ • f) - Sq • Sq/3 
with f = r/r, S = Sq -|- Sg and L = r x (-iV) ptT] . 
Here, V , Vs, Vt and Vls represent the spin- independent 
central, spin-spin, tensor and spin-orbit potentials, re- 
spectively. 

In this paper, we focus only on the S'-wave charmo- 
nium states (77c and J/tp). We perform an appropriate 
projection with respect to the discrete rotation, which 
provides the BS wave function projected in the A^ rep- 
resentation, (/>r(i') — ^ cj)T{A'^\r). This projected BS wave 
function corresponds to the S'-wave in continuum theory 
at low energy [TH]. We simply denote the A'^ projected 
BS wave function by (f>rir) hereafter. 

The stationary Schrodinger equation for the projected 
BS wave function (prii") is reduced to 



+ V{r) + Sq ■ S^Vsir) Mr) = i?r0r(r) (4) 



at the leading order of the v-expansion [TS] . The spin op- 
erator Sq • Sg can be easily replaced by expectation val- 
ues —3/4 and 1/4 for the PS and V channels, respectively. 
As a result, both spin-independent and -dependent part 
of the central interquark potentials can be separately 
evaluated through a linear combination of Eq. (|4| calcu- 
lated for both PS and V channels as 



V{r) 
Vs{r) 



1 f3 V^0v(r) 
mq I4 
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'hyp 



Mr) ' 4 Mir) f^^' 
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mg \ <j)vir) Mir) 



(6) 



where Save = A/avc ^ 2toq and £^hyp = A/y — Mps. The 
mass Mavo denotes the spin- averaged mass as ^Mps -I- 
|Mv. The derivative is defined by the discrete Lapla- 
cian with nearest-neighbor points. As for other spin- 
dependent potentials such as the tensor potential Vt and 
the spin-orbit potential VlSj in principle, this approach 
can access them by considering the P-wave quarkonium 
states such as the Xc (0^^, 1''"^) and he (1^ ) states, 
which must leave contributions of Vt and Vls to Eq.(|4|. 

It is worth mentioning that the quark kinetic mass 
mg, which is essentially involved in the definition of the 
interquark potentials, can be self-consistently evaluated 
within the same framework [3]. The proper determina- 
tion of the quark mass has a key role in establishing the 
connection to the static heavy-quark potential given by 



Wilson loops in the infinitely heavy-quark limit. A more 
detailed discussion can be found in Ref. [9]. 

To calculate the charmonium potential, we have per- 
formed dynamical lattice QCD simulations on a lattice 
L3 X T = 323 X 54 ^j^j^ 2+1 flavor PACS-CS gauge config- 
urations generated by Iwasaki gauge action at /3 = 1.90, 
which corresponds to a lattice cutoff of « 2.2 GeV 
(a « 0.091fm) [1^. The spatial lattice size corresponds 
to L « 3 fm. The hopping parameters for the light 
sea quarks {Kud,«^s}={0. 13781, 0.13640} correspond to 
= 156(7) MeV and Mk = 554(2) MeV [12]. Al- 
though the light sea quark masses are slightly off the 
physical point, the systematic uncertainty due to this fact 
could be extremely small in our project. Our results are 
analyzed on all 198 gauge configurations, which are avail- 
able through International Lattice Data Grid and Japan 
Lattice Data Grid [20 . We fix gauge configurations to 
Coulomb gauge. 

For the charm quark, we employ the RHQ action to 
control systematic uncertainties coming from the dis- 
cretization error induced by large quark mass [13]. The 
RHQ action utilized here is a variant of the Fermilab ap- 
proach [21j and has five parameters Kc, i^, ^g, cb and ce- 
The parameters Ts, cb and ce are determined by tad- 
pole improved one-loop perturbation theory |22j . For 
we use a nonperturbatively determined value, which is 
adjusted by reproducing the effective speed of light to be 
unity in the dispersion relation £^^(p^) = + c^gjpp 
for the spin-averaged IS* charmonium state, since the 
parameter v is sensitive to the size of hyperfine mass 
splitting [23]. We choose Kc to reproduce the exper- 
imental spin-averaged mass of IS" charmonium states 
Mf^P = 3.0678(3) GeV. To calibrate adequate RHQ pa- 
rameters, we employ a gauge invariant Gauss smearing 
source for the standard two-point correlation function 
with four different finite momenta. As a result, the rel- 
evant speed of light, c^g = 1.04(5), is observed for the 
spin-averaged mass of IS' charmonium states with our 
chosen RHQ parameters summarized in Table [l| 

We have computed charm quark propagators with two 
wall sources located at different time slices tg/a = 6 and 
57, to increase statistics. Dirichlet boundary conditions 
are imposed for the time direction to eliminate unwanted 
contributions across time boundaries. We calculate a pair 
of four-point correlation functions from two wall-source 
quark propagators and fold them together to create a 
single four-point correlation function. 

Low-lying S- and P-wave charmonium masses mea- 
sured in this study are all close to experimental values, 
though the hyperfine mass splitting Mhyp = 0.1124(9) 
GeV is slightly smaller than the experimental value, 
M^yP = 0.1166(12) GeV. Note that we simply neglect 
the disconnected diagrams in both the rjc and J/ip cor- 
relation functions. The similar value of the hyperfine 
mass splitting is reported even on the physical point in 
Ref. [23;. We summarize resulting charmonium masses 
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FIG. 1: Spin-independent interquark potentials calculated 
from both the BS wave function of the charmonium states (+ 
and •) and the standard Polyakov line correlator (x and ■). 
Filled symbols represent on-axis data. The constant terms 
are subtracted from both potentials to set to V{ro) — 0. 

in Table H 

TABLE II: Summary of the Cornell parameters and the quark 
mass determined from lattice QCD. For comparison, the cor- 
responding values adopted in a NRp model [5] are also in- 
cluded. 

This work Polyakov lines NRp model 
on-axis full set 
A 0.861(17) 0.813(22) 0.403(24) 0?7281 
^ [GeV] 0.394(7) 0.394(7) 0.462(4) 0.3775 
niQ [GeV] 1.74(3) oo 1.4794 



First, we show a result of the spin-independent charmo- 
nium potential V{r) in Fig. [ij where the constant term 
is subtracted to set V{ro) = with the Sommer scale. 
To ~ 0.5 fni. For comparison, the static heavy-quark 
potential calculated from the Polyakov line correlator, 
which corresponds to the one obtained in the infinitely 
heavy-quark limit similar to the Wilson loop results [5], 
is also displayed in Fig. [ij As expected, the charmo- 
nium potential calculated in the BS amplitude method 
properly exhibits the linearly rising potential at large dis- 
tances and the Goulomb-like potential at short distances. 

Here we give some technical remarks on systematic 
uncertainties. In the BS amplitude method, we take a 
weighted average of data points in the wide range of 
\t — ts\/a = 26 - 48 for determining the equal-time BS 
wave function. Therefore, the resulting charmonium po- 
tential has a much smaller systematic error stemming 
from the uncertainty in the choice of time window than 
the conventional approach to calculate the static heavy- 
quark potential by Wilson-loops or Polyakov lines. On 
the other hand, the discretization error seems to much 
severely appear in the charmonium potential especially 
near the origin. The Goulomb-like behavior obtained in 
the BS amplitude method may contain large uncertain- 
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FIG. 2: Spin-independent charmonium potential calculated 
from the BS wave function. The dashed curve is the fit- 
ting result by the Cornell parametrization. The shaded band 
shows statistical fitting uncertainty calculated by the jack- 
knife method. For comparison, the phenomenological poten- 
tial adopted in a NRp model [5] is also included as solid curve. 
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FIG. 3: Spin-spin charmonium potential calculated from 
the BS wave function. The dashed, dotted and dash-dotted 
curve correspond to fitting results of Yukawa, exponential and 
Gaussian functional forms, respectively. For comparison, the 
phenomenological potential adopted in a NRp model [5] is 
also included as solid curve. 



ties, which should vanish in the continuum hmit. To 
avoid the large discretization error, we hereafter prefer 
to use the "on-axis" data, which less suffers from the 
rotational symmetry breaking in the finite cubic box. 

The Cornell parametrization is simply adopted for fit- 
ting our data of the spin-independent central potential: 

V{r) = -- +ar + Vo (7) 
r 

with the Coulombic coefficient A, the string tension a 
and a constant Vq- We have carried out correlated fits 
with full covariance matrix for on-axis data over range 
4 < r/a < 10, while uncorrelated fits are adopted in 
full data analysis including all off-axis data points due 
to high correlation between different r points. The fit- 
ting results are listed in Table |l] together with the phe- 
nomenological values employed by a nonrelativistic po- 
tential (NRp) model in Ref. [S]. From on-axis data only, 
we get the Cornell parameters of the charmonium po- 
tential: A = 0.861(17) and = 0.394(7) MeV with 
acceptable x^/dof (« 2.2). The quoted errors represent 
only the statistical errors given by the jack-knife analysis. 

In Fig. [2] we show on-axis data points of the spin- 
independent charmonium potential with the fitted curve 
(dashed curve). The phenomenological potential used in 
NRp models [5] is also plotted as a solid curve for com- 
parison. As shown in Fig. [2j although the charmonium 
potential obtained from lattice QCD is quite similar to 
the one in the NRp models, the string tension of the 
charmonium potential is slightly stronger than the phe- 
nomenological one. Therefore our result indicates that 
there are only minor modifications required for the spin- 
independent central potential in the NRp models. 

Moreover, it seems that a gap for the Coulombic co- 
efficients between the conventional static potential from 



Wilson-loops and the phenomenological potential used in 
the NRp models closes by our new approach, which non- 
perturbatively accounts for a finite quark mass effect. 

It is worth mentioning that the string breaking, which 
would be induced by the presence of dynamical quarks 
was not observed at least in the range r < 1 fm, where we 
still get a better signal-to-noise ratio. It is indeed what 
we expected, since we cannot access information of the 
potential outside of the localized wave function, which 
represents the charmonium bound state within the BS 
amplitude method. We here calculate only the BS wave 
functions of 15 charmonium states, which are quickly 
dumped around outside of r > 1 fm. Therefore, at least 
the similar calculation for the higher-lying charmonium 
states, whose wave function can be extended until the 
string breaking sets in, is demanded to observe such ef- 
fect. 

In this calculation, the kinetic mass of the charm quark 
is determined self-consistently within the BS amplitude 
method as well. (See Ref. [9 for details.) The charm 
quark mass obtained in this study is about 17% heavier 
than the one adopted in the NRp models, of which value 
is also listed in Table HIl This difference should not be 
taken seriously since the spatial profile of the spin-spin 
potential from lattice QCD is slightly different from the 
one used in the NRp models as we will discuss later. 

In Fig. [Sj we show the spin-spin term of the charmo- 
nium potential and the corresponding phenomenological 
one found in Ref. Our spin-spin potential exhibits 
the short-range repulsive interaction, which is required 
by the charmonium spectroscopy, where the higher spin 
state in hyperfine multiplets receives heavier mass. It 
should be reminded that the Wilson loop approach fails 
to reproduce the correct behavior of the spin-spin interac- 
tion, since the leading-order spin-spin potential classified 
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TABLE III: Results of fitted parameters for the spin-spin po- 



tential with three 


types of the fittin 


g functional form 




functional form 


a 


P 


xVdof 


Yukawa-type 


0.287(8) 


0.894(32) GeV 


7.28 


Exponential- type 


0.825(19) GeV 


1.982(24) GeV 


1.46 


Gaussian-type 


0.314(4) GeV 


1.020(11) GeV^ 


22.79 



in potential nonrelativistic QCD becomes attractive at 
short distances [8 . 

In contrast of the case of the spin-independent poten- 
tial, the spin-spin potential obtained here is absolutely 
different from a repulsive (S-function potential generated 
by perturbative one-gluon exchange, which is widely 
adopted in the NRp models. However, such the con- 
tact form c>c i5(r) of the Fermi-Breit type potential is not 
reliable since the pointlike spin-spin interaction can not 
give a finite hyperfine mass splitting of the P- or higher- 
wave charmonia p^. Indeed, the finite-range spin-spin 
potential described by the Gaussian form is adopted in 
Ref. [51 , where many properties of conventional charmo- 
nium states at higher masses are predicted. 

This phenomenological spin-spin potential is also plot- 
ted in Fig. [3] for comparison. Although there is a slight 
difference at very short distances, where systematic un- 
certainties become severe due to lattice discretization er- 
rors, the spin-spin potential from first principles QCD is 
barely consistent with the phenomenological one. 

To examine the appropriate functional form for the 
spin-spin potential, we have tried three types of func- 
tional forms: 



Vs{r) 



a exp(— /3r)/r 
a exp(— /3r) 
a exp(— /3r^) 



Yukawa form 
Exponential form 
Gaussian form. 



(8) 



We then determine which functional form can give a rea- 
sonable fit over the range of r/a from 2 to 10. All results 
of correlated fits are summarized in Table III The 



long-range screening observed in the spin-spin potential 
is more easily accommodated by the Yukawa form or the 
exponential form than the Gaussian form that is often 
employed in the NRp models. Although the exponential 
form provides the smaller ;\;^/dof than the Yukawa form, 
a solid conclusion requires more accurate information on 
the short-range behavior of the spin-spin potential. 

In this paper, we have studied both spin-independent 
and spin-dependent part of the charmonium potential 
by means of the BS wave function of 15 charmonium 
states in dynamical lattice QCD simulations. The spin- 
independent charmonium potential obtained from lattice 
QCD with almost physical quark masses is quite simi- 
lar to the one used in the NRp models. The spin-spin 
potential, which is, for the first time, determined in dy- 
namical lattice simulations, exhibits not pointlike, but 
finite-range repulsive interaction. Its r-dependence is 



barely consistent with the phenomenological one adopted 
in Ref. [5j. Thus, a full set of the reliable spin-dependent 
potentials derived from lattice QCD within our approach 
can provide new and valuable information to the NRp 
models. This improvement of the spin-dependent po- 
tentials will help in making accurate theoretical predic- 
tions about the higher-mass charmonium states. We plan 
to extend our research to determine all spin-dependent 
terms in the charmonium potential, including the tensor 
and spin-orbit forces and also to address all the possible 
systematic uncertainties described in the text. 
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